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Abstract. We introduce a new domain for finding precise numerical invariants of pro- 
grams by abstract interpretation. This domain, which consists of sub-level sets of non- 
linear functions, generalizes the domain of linear templates introduced by Manna, Sankara- 
narayanan, and Sipma. In the case of quadratic templates, we use Shor's semi-definite 
relaxation to derive safe and computable abstractions of semantic functionals, and we 
show that the abstract fixpoint can be over-approximated by coupling policy iteration and 
semi-definite programming. We demonstrate the interest of our approach on a series of 
examples (filters, integration schemes) including a degenerate one (symplectic scheme). 



We introduce a complete lattice consisting of sub-level sets of (possibly non-convex) func- 
tions, which we use as an abstract domain in the sense of abstract interpretation |CC77j for 
computing numerical program invariants. This abstract domain is parametrized by a basis 
of functions, akin to the approach put forward by Manna, Sankaranarayanan, and Sipma 
(the linear template abstract domain |SSM05j see also [SCSM06J), except that the basis 
functions or templates which we use here need not be linear. The domains obtained in this 
way encompass the classical abstract domains of intervals, octagons and linear templates. 
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X 


= [0,1]; 




V 


= [0,1]; 




h 


= 0.01; 




w 


hile (true) 


{ [2] 




w = v ; 






v = v*(l — h) 


— h*x ; 




x = x+h*w; 


[3] } 




{x 2 < 3.5000, v 2 < 2.3333, 2x 2 + 3w 2 + 2xv < 7} 



Figure 1: Euler integration scheme of a harmonic oscillator and the loop invariant found at 
control point 2 



To illustrate the interest of this generalization, let us consider a harmonic oscillator: 
x + cx + x = 0. By taking an explicit Euler scheme, and for c = 1 we get the program 
shown at the left of Figure [TJ 

The invariant found with our method is shown right of Figure [T] For this, we have 
considered the template based on functions {x 2 ,v 2 ,2x 2 + 3v 2 + 2xv}, i.e. we consider a 
domain where we are looking for upper bounds of these quantities. This means that we 
consider the non-linear quadratic homogeneous templates based on {x 2 , v 2 }, i.e. symmetric 
intervals for each variable of the program, together with the non-linear template 2x 2 + 
3v 2 + 2xv. The last template comes from the Lyapunov function that the designer of 
the algorithm may have considered to prove the stability of his scheme, before it has been 
implemented. This allows us to represent set defined by constraints of the form x 2 < c±, 
v 2 < C2 and 2x 2 + 3v 2 + 2xv < C3 where ci, C2 and C3 are degrees of freedom. In view 
of proving the implementation correct, one is naturally led to considering such template^] 
Last but not least, it is to be noted that the loop invariant using intervals, zones, octagons 
or even polyhedra (hence with any linear template) is the insufficiently strong invariant 
h = 0.01 (the variables v and x cannot be bounded.) However, the main interest of the 
present method is to carry over to the non-linear setting. For instance, we include in our 
benchmarks a computation of invariants (of the same quality) for an implementation of a 
highly degenerate example (a symplectic integration scheme, for which alternative methods 
fail due to the absence of stability margin). 

Contributions of the paper. We describe the lattice theoretical operations in terms of Galois 
connections and generalized convexity in Section [2] We also show that in the case of a 
basis of quadratic functions, a good over-approximation F^~ of the abstraction F$ of a 
semantic functional F can be computed by solving a semi-definite program (Section [3]). 
This over-approximation is obtained using Shor's relaxation. The advantage of the latter 
is that it can be solved in polynomial time to an arbitrary prescribed precision by the 
ellipsoid method [GLS88J , or by interior point methods |NN94j if a strictly feasible solution 
is available (however, we warn the reader that interior points methods, which are more 
efficient in practice, are known to be polynomial only in the real number model, not in the 
bit model, see the survey |PR97j for more information). 



^Of course, as for the linear templates of SSM05, SCSM06 , we can be interested in automatically finding 
or refining the set of templates considered to achieve a good precision of the abstract analysis, but this is 
outside the scope of this article. 
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Moreover, we show in Subsection 4.3 that the vectors of Lagrange multipliers produced 
by this relaxation correspond to policies, in a policy iteration technique for finding fixpoints 
or at least postfixpoints of F^, precisely over-approximating the fixpoints of FK Finally, we 
illustrate on examples (linear recursive filters, numerical integration schemes) that policy 
iteration on such quadratic templates is efficient and precise in practice, compared with 
Kleene iteration with widenings/narrowings. The fact that quadratic templates are efficient 
on such algorithms is generally due to the existence of (quadratic) Lyapunov functions that 
prove their stability. The method has been implemented as a set of Matlab programs. 

Related work. This work is to be considered as a generalization of [SSM05], [SC SM06] 
because it extends the notion of templates to non- linear functions, and of |CGG + 05] , 
|GGTZ07| . |AGG08| . |GS07a| and |GS07bj since it also generalizes the use of policy itera- 
tion for better and faster resolution of abstract semantic equations. Polynomial inequalities 
(of bounded degree) were used in [BRCZ05] in the abstract interpretation framework but 
the method relies on a reduction to linear inequalities (the polyhedra domain), hence is 
more abstract than our domain. Particular quadratic inequalities (involving two variables - 
i.e. ellipsoidal methods) were used for order 2 linear recursive filters invariant generation in 
[Fcr05pl Polynomial equalities (and not general functional inequalities as we consider here) 
were considered in [MOS04, RCK07J. The use of optimization techniques and relaxations 
for program analysis has also been proposed in |Cou05| . mostly for synthesising variants 
for proving termination, but invariant synthesis was also briefly sketched, with different 
methods than ours (concerning the abstract semantics and the fixpoint algorithm). Finally, 
the interest of using quadratic invariants and in particular Lyapunov functions for proving 
control programs correct (mostly in the context of Hoare-like program proofs) has also been 
advocated recently by E. Feron et al. in [FF081 lFA"U8] . 

Finally, we note that a preliminary account of the present work appeared in the con- 
ference paper [AGG10] . 

2. Lattices of sub-level sets and P-support functions 

We introduce a new abstract domain, parametrized by a basis of functions (P below). The 
idea is that an abstract value will give the bounds for each of these functions, hence the 



name of sub-level sets, with some abstract convexity condition, Definition 2.9 The abstract 
values are computed as the supremum of each function of the basis on a certain sub-level 
set. The name support functions is due to the similarity with the classical support function 
from convex analysis where the functions of the basis are all linear. The idea of using the 
classical support function in system verification appeared independently in [LGG09] . 

2.1. P-sub-level sets, and their Galois connection with V(M. d ). Let P denote a set 
of functions from R d to M, which is going to be the basis of our templates. The set P is 
not necessarily finite and the functions p £ P are not necessarily linear. We denote by 
T(P,R) tn e set of functions v from P to 1 = R U {±oo}. We equip P(P,I) with the 
classical partial order for functions i.e. v < w <J=^ v(p) < w(p) for all p G P. We 
order the set of a ll subsets of R d by the subset relation C. We define a Galois connection 



(Proposition 2.5) between J-(P, K) and the set of subsets of R d (made of a concretisation 



2 A generalization to order n linear recursive filters is also sketched in this article. 
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Figure 2: A P-sub-level set aris- 
ing from non-convex qua- 
dratic functions. 



{!/-«< 3, y + x<3, -:y<0} 



Figure 3: A P-sub-level set arising 
from linear forms. 



operator v i— > v*, Definition 2.1 and an abstraction operator C i— > C', Definition 2.4). This 



will give the formal background for constructing abstract semantics using P-sub-level sets 
using abstract interpretation |CC77| . in Section [3j 

Definition 2.1 (P-sub-level sets). To a function v £ P(P, M), we associate the P-sub-level 
set denoted by v* and defined as: 

v* = {x £ R d | p(x) < v(p), Vp G P} 

The notion of sub- level set is well known in convex analysis. When P is a set of function 
that are convex, the P-sub-level sets are convex. In our case, P can contain non-convex 
functions so P-sub-level sets are not necessarily convex in the classical sense. 

Example 2.2. We come back to the first example and we are focusing on its representation in 
terms of P-sub-level set. Let us write, for (x, v) £ M 2 , p\ : (x, v) \-t x 2 , p2 ■ (x, v) \-t v 2 and 
p 3 : (x,v) i-> 2x 2 + 3v 2 + 2xv. Let us take P = {pi,p 2 ,p 3 }, w(pi) = 3.5000, v(p 2 ) = 2.3333 
and v(ps) = 7. The set v* is precisely the one shown right of Figure [TJ 

Example 2.3. We next show some P-sub-level sets which are not convex in the usual sense. 
Let us write, for (x, y) £ M 2 , p\ : (x, y) i-> — y 2 — (x + 2) 2 , p 2 : (x, y) \-¥ —y 2 — (x — 2) 2 and 
p 3 : (x,y) (->■ -(y - 2) 2 - x 2 , p 4 : (x,y) (->■ -(y + 2) 2 - x 2 . Let us take P = {pi,P2,P3,P4} 
and v(pi) = v(p 2 ) = v(ps) = v(p^) = —2. The set v* is shown Figure [2] 

In our case, P is a set of functions from M. d to R possibly non-linear, so we generalize 
the concept of support functions (e.g. see Section 13 of |Roc96j ). The idea of generalizing 
support functions to the non-linear case is not new and this extension is due to Moreau 
|Mor70j . 

Definition 2.4 (P-support functions). To X C 1^, we associate the P-support function 
denoted by X^ and defined as: 

X'(p) = sup^(x) 

Proposition 2.5. The pair of maps v i— >■ v* and X t— > X* defines a Galois connection 
between P(P, M) and the set of subsets ofM. d . 

Proof. First, the functions and X i— >■ X* are clearly monotonic. Now, we have to 

show that, for all X C R d and all v £ P(P,I): X C v* <^> Xt(p) < «(p), Vp G P. 
Let X Q W l , v £ P(P,1) and p £ P. We assume that X C «*. This implies that 
X'(p) = sup{p(x) | x G X} < sup{p(x) | q(x) < v(q), \/q £ P}. But, for every x G v*, 
p(x) < v(p), so sup{p(x) I q(x) < v(q), \/q £ P} < v(p), hence, X*(p) < v(p) and the first 
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{y 2 - x 2 < 9, x < 3, -x < 3} 



{x — y < 3, y + x < 3, —x < 3} 

Figure 5: p2-convex hull of example 

Figure 4: Pi-convex hull of example 



implication is shown. Now, we assume that X^(p) < v(p), Vp G P. Let i £ X, we have, 
for all p G P, p(x) < X^(p), thus p(x) < v(p) and finally iB*. □ 

Remark 2.6. In the proof, we showed < t>, we can similarly prove that X C (X')*. 

Using this remark and using the monotonicity of v i— )■ u* and X i— > X', we get u* C 
((w*)')* C and similarly, {(X')*y = X*. These properties are well-known in Galois 
connection theory. 

In the terminology of abstract interpretation, is the abstraction function, usually 
denoted by a, and (.)* is the concretisation function, usually denoted by 7. 



2.2. The lattices of P-convex sets and P-convex functions. The sets of points in M d 
which are exactly represented by their corresponding P-sub-level sets are called P-convex 
sets, as in the definition below. These can be identified with the set of abstract elements 



we are considering^! We show in Theorem 2.12 that they constitute a complete lattice. We 



use the terminology P-convex because of the analogy with the abstract convexity defined 
by Moreau |Mor70j . Singer |Sin97| and Rubinov [RubOOj. These authors define convexity 
without linearity from a given family of (non-linear) functions. Moreau called this notion 
of convexity, P-convexity where B is the family of functions whereas Singer and Rubinov 
use the term abstract convexity. The P-convexity defined in this paper corresponds to 
the classical notions of closure from the Galois connections theory and also to the abstract 
convexity introduced by Moreau. 

Definition 2.7 (P-convex hull). The P-convex hull of an element v G P(P, M) is the 
function vexp(w) which is equal to (v*y. Similarly, the P-convex hull of a subset X is the 
set vexp(X) which is equal to (X^)*. 

Example 2.8. Let us consider the triangle, depicted in Figure[3} Let us take Pi = {(x, y) i-> 
y + x, (x, y) i->- x — y, (x, y) h-> — x}. Its Pi-convex hull is the one depicted Figure |4j If we 
take instead P2 = {(x, y) y 2 — x 2 , (x, y) h-> x, (x, y) :h-» — x}, its P2-convex hull is shown 
in Figure [5j 

^Formally, this is the upper-closure in "P(R d ) of the set of abstract elements. 
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Definition 2.9 (P-convexity). Let v S P(P, M), we say that v is a P-convex function if 
v = vexp(u ). A set X C W 1 is a P-convex set if X = vexp(X). 

Example 2.10. Let us come back to the triangle depicted in Figure [3] If P is the set of 
linear forms defined by the faces of this triangle i.e. P consists of the maps (x, y) :t— >■ y — x, 
(x, y) :i— > y + x and (x, y) :i— )■ — y, then it is a P-convex set. But if P is, for example, linear 
forms defined by orthogonal vectors to the faces of the triangle, the previous triangle is no 
longer an P-convex set. 

Abstract (P-) convexity is a special instance of Galois connection. The latter are classi- 
cally used in abstract interpretation. In particular, the theory of Galois connections yields 
the following result: the P-convex hull of a function v G P(P, M) is the greatest P-convex 
function which is smaller than v and dually the P-convex hull of a subset X is the smallest 
P-convex set which is greater than X. 

We respectively denote by Vexp(P i— >■ M.) and Vexp(M <i ) the set of all P-convex function 
of P(P, M) and the set of all P-convex sets of M. d . 

Definition 2.11 (The meet and join). Let v and w be in P(P, K). We denote by inf(v,w) 
and sup(f , w) the functions defined respectively by: 

p i->- mi(v(p),w{p)) and p h-» sup(v (p),w(p)) . 
We equip Vexp(P i-> M) with the meet (respectively join) operator: 

v V w = swp(v, w) (2-1) 
v A u; = (inf(u,u;)*) t (2.2) 
Similarly, we equip Vexp(M d ) with the two following operators: X U Y = ((X U Yy)*, 

xnY = xnY. 

The family of functions Vexp(P i— > M) is ordered by the partial order of real- valued 
functions i.e. v < w <^=^ v(p) < w(p) Vp € P. The family of set Vexp(JR <J ) is ordered by 
the subset relation denoted by C. The next theorem follows readily from the fact that the 
pair of functions v i— >■ v* and X i— )■ defines a Galois connection, see e.g. [DP021 § 7.27]. 

Theorem 2.12. (Vexp(P i-> K),A, V) and (Vexp(]R <1 ), n, U) are isomorphic complete lat- 
tices, n 



2.3. Intervals, Zones, Octagons and Sankaranarayanan et al.'s linear templates. 

The domain of intervals arises as a special domain of P-convex sets, in which the basis is 
P = {x%, —x\, . . . , x n , —x n } where Xi (i = 1, . . . , n) are the program variables. An abstract 
value v in our domain encodes the supremum v{xj) and the infimum — v(— Xi) of an interval 
for the variable X{. 

Zones and octagons are treated in a similar manner. For instance, for zones, take 
P = {Xi — Xj | i,j — 0, . . . , n, i 7^ j}, adding a dummy variable xq (always equal to 0), 
as customary, to subsume intervals. Of course, linear templates as defined in |SSM05| are 
particular P-convex sets, for which P is given by a finite set of linear functionals. 

We remark that in the case of zones, v(x{ — Xj) is exactly the entry i, j of the DBM (Dif- 
ference Bound Matrix) representing the corresponding zone. Also, elements of Vexp(P i— > M) 
corresponding naturally to closed DBMs, that is, canonical forms of DBMs. As is well known 
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[Min(2], the union of two zones preserves closure whereas the intersection does not neces- 
sarily preserve closure. More generally, the set of closed elements for a Galois connection is 
stable under the join operation. This is reflected in our domain by (2.1) and (2.2). 



3. Quadratic templates 

In this section, we instantiate the set P to linear and quadratic functions. This allows us to 
give a systematic way to derive the abstract semantics of a program. The main result is that 
the abstract semantics for an assignment and for a test, can be safely over-approximated 
by Shor's relaxation scheme, Theorem |3.13| 

Definition 3.1. We say that P is a quadratic zone iff every element template p £ P can 
be written as: 

x i Y p(x) = x T A p x + bpX, 

where A p is a d x d symmetric matrix (in particular A p can be a zero matrix), x T denotes 
the transpose of a vector x, b p is a M. d vector. 

Now, we suppose that P is finite. We denote by F(P, M) the set of functions from P 
to E, P(P,RU {+oo}) the set of functions from FtolU {+00} and F(P,R + ) the set of 
functions from P to M.+ . 



Suppose now we are given a program with d variables (x\, . . . , Xd) and n control points 
numbered from 1 to n. We suppose this program is written in a simple toy version of 
a C-like imperative language, comprising global variables, no procedures, assignments of 
variables using only parallel assignments {x\, . . . ,Xd) = T(x±, . . . , Xd), tests of the form 
(xi, . . . , Xd) £ C, where C is some shape in V(M. d ), and while loops with similar entry tests. 
We do not recapitulate the standard collecting semantics that associates to this program a 
monotone map F : (V(lR. d )) n — > (V(M. d )) n whose least fixpoints ffp(-F) has as ith component 
(i = 1, . . . , n) the subset of R d of values that the d variables x%,...,Xd can take at control 
point i. 

The aim of this section is to compute, inductively on the syntax, the abstraction (or a 
good over-approximation of it) of F from P(P, M) n to itself defined as usual as: 

F\v) = (F(v*Y) (3.1) 

The notation v* is in fact the vector of sets (v*, ■ ■ ■ , u*) and (F(v*y) is also interpreted 
component- wise. The notation vexp(f) will be also understood component- wise. 



3.1. Shor's semi-definite relaxation scheme. Finding the maximal value of a non- 
concave quadratic function under convex or non-convex quadratic constraints is known to 
be an NP-Hard problem, see [Vav90] for a discussion of complexity issues in quadratic 
programming. Shor's relaxation scheme (see [TNOl} Section 4.3.1] or Shor's original ar- 
ticle |Sho87| for details) consists of over-approximating the value of a general quadratic 
optimization problem by the optimal value of a semi-definite programming (SDP for short) 
problem, the latter being computationally more tractable. 

Indeed, SDP problems can be solved in polynomial time to an arbitrary prescribed 
precision by the ellipsoid method [GLS88J. More precisely, let e > be a given rational, 
suppose that the input data of a semi-definite program are rational and suppose that an 
integer ./V is known, such that the feasible set lies inside the ball of the radius ./V around 
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zero. Then an e-optimal solution (i.e., a feasible solution the value of which is at most 
at a distance e from the optimal value) can be found in a time that is polynomial in the 
number of bits of the input data and in — log e. Moreover, an e-solution of an SDP problem 
can also be found in polynomial time by interior point methods [NN94J if a strictly feasible 
solution is available. However, when the input is rational, no size on the bit lengths of the 
intermediate data is currently known, so that the term "polynomial time" for interior point 
methods is only understood in the model of computation over real numbers (rather than 
in bit model [GJ79]). The advantage of interior methods is that they are very efficient in 
practice. We refer the reader to [P R97j for more information. 

Let /, {/i}i=i,...,m be quadratic functions on M. d . Let us consider the following con- 
strained maximization problem: 

sup{/(x) I fi(x) <0, Vi = l,...,m} (3.2) 

In constrained optimization, it is classical to construct another constrained optimization 
problem from the initial one in order to solve an easier problem. A technique called Lagrange 
duality (for details see for example [AT031 Section 5.3]) consists in adding to the objective 
function the inner product of the vector of constraints with a positive vector of the euclidean 



space whose the dimension is the number of constraints. In our context, the value of (3.2) 
is given by the following sup-inf (primal) value (|3.3[): 



sup A f m ~ X ifi( x ) ■ ( 3 - 3 ) 
A simple result of constrained optimization called weak duality theorem ensures that 



if we commute the inf and the sup in the formula (3.3), the result is greater than the value 



(3.3). The commutation of the inf and the sup gives us the so called (dual) value: 

m 

inf sup f{x)- hfi{x) ■ (3.4) 



The vectors A £ are called vectors of Lagrange multipliers. The function A i— > 
su PxeiR d f( x ) ~ S2=i ^ifi( x ) is always convex and lower semi-continuous (these definitions 
are recalled in the appendix), so it has good properties to minimize it. If the function / 
is concave, the functions fi are convex and if the Slater constraint qualification (i.e. there 



exists x 6 Mr such that fi(x) < for all i = 1, . . . ,m) holds then (3.3) and (3.4) coincide 



(this will be used Proposition 3.8). 



Shor's relaxation scheme consists in computing the value (3.4) by solving a semi-definite 
program. We introduce the matrix M(g), for a quadratic function g written as x T A g x + 
bjx + c g and the matrix N(y) for a real y defined as: 

M{g) = ^ , and N ljX {y) = y, iV M (y) = if (i,j) ? (1, 1) (3.5) 

Let -< denote the Lowner ordering of symmetric matrices, so that A -4 B iff all eigen- 
values of B — A are non-negative. 

When we fix A € W+, we have to maximize an unconstrained quadratic problem and the 
the maximum is finite iff there exists r\ G M such that, for all x £ fix)— ~Y^h=\ ^ifi( x ) — V 
and since /, fi are quadratic functions, this is equivalent to x T {Af — Y^iLi ^i-^fi) x + (&/ — 
YliLi Kbfc + Cf — YliLi ^i c fi — < for all x G ~R d which is equivalent to the fact that the 
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matrix M(f) + r)N(—l) — YliLi \M(fi) is negative semi-definite. Consequently, taking the 



infimum over A G WV:, we recover the value (3.4). 



In conclusion, Shor's relaxation scheme consists in solving the following SDP problem: 

rn 

Mm 77 s.t. M{f)+r 1 N{-l)-Y J \M{f i )]<Q (3.6) 



i=l 



which is equal to the value (3.4), hence an over-approximation of the optimal value of the 



problem (3.2). We can use a verified SDP solver as VSDP [JCK07J to solve a SDP problem. 



Example 3.2. We consider the rotation of angle (ft G (0, 7r) leaving invariant the unit circle 
S 1 = {(x,y) £ R 2 \ x 2 + y 2 = 1}. We write T the linear map associated to the rotation 
matrix A: 

'A =A{x , y f = H<t> ~ sin /]( x 

y J ysm <p cos <p J \y 

where x 2 + y 2 = 1 . To show that the proposed relaxation is accurate we next observe that 
it preserves the fact that the unit sphere is invariant by rotation i.e. T(S 1 ) = S . The 
inclusions S 1 C T(S 1 ) and T(S' 1 ) C S 1 can be proved by the same manner because A is 
an orthogonal matrix i.e. A T A = AA T = Id: ((x,y) G 5 1 ==> (x,y) G T(S 1 )) <J=^ 
(A T (x, y) G S 1 for {x, y) G S 1 ), hence, to show S 1 C T(5 : ) is reduced to show T^S 1 ) C S 1 
where T* is the linear map associated to ^4 T . We only prove that T(S' 1 ) C S* 1 and we use 
the Shor's relaxation scheme to prove that. 

We introduce the set of quadratic functions P = {pi(x, y) H > x 2 + y 2 , p 2 (x,y) h-> 
— (x 2 + y 2 )} and we set i>i(pi) = 1 and ^1(^2) = — 1- The unit sphere is the P-sub- level set 
of v%\ {{x,y) G R 2 I p(x,y) < v\(p), Vp G P}. To show T(5 X ) C S , it suffices to show 
T"(u*)t = wi. We write: 

«2(pi) := ^WJ^Pi) = sup{pi (T(x,y)) \ pi(x,y) < 1, p 2 (x,y) < -1} 
^2(^2) := 3 n (v^) t (p 2 ) = sup{p 2 (T(x,y)) | pi{x,y) < 1, p 2 (x,y) < -1} 



For the calculus of v 2 (pi), the matrices defined in 3.5 



arc: 



M(pioT)= 1 , M(p 2 oT) = -M(pioT), 

\0 1/ 

/o o\ 

M(pi) = 1 and M(p 2 ) = -M(pi) . 



By using the Shor's relaxation scheme, the SDP problem (3.6) gives the following equalities: 



v 2 {pi) = Min 77 s.t. M(p 1 oT)+r ] N(-l)+X(p 1 )(N(l)-M(p 1 ))+X(p2)(N(-l)-M(p2)) ^ 

A(pi)>0 
A(p 2 )>0 

and 

v 2 (p 2 ) = Min 77 s.t. M(p 2 oT)+r/iV(-l)+A(pi)(^(l)-M(pi))+A(p2)(A r (-l)-Af(p2)) =< 

A(pi)>0 
A(p 2 )>0 
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We can rewrite the two previous SDP problems as follows: 

f -V + Hpi) ~ HP2) o 

v 2 (pi)= Min rys.t. 1 - A(p x ) + A(p 2 ) |^0 



and 

{ -v + Kpi) - Hp*) o 

w 2 (p 2 ) = Min r/s.t. -1 - X( Pl ) + X(p 2 ) |^0 

ZMl V -1-AfeO + AfoO 

To solve these optimization problems, we could call an SDP solver, but in this case, it 
suffices to solve a system of inequalities: all the diagonal elements must be non-positive, for 
example, for the first problem, this implies that r/ > 1 and since we minimize 77 we get 77 = 1. 
Hence, we find v^ipi) = 1 and ^2(^2) = —1 and finally we have proved = v± and we 

conclude that the rotation invariance of the unit sphere S* 1 = {(x, y) 6 M? \ x 2 + y 2 = 1} is 
preserved for the relaxation. 



3.2. Abstraction of assignments and tests using Shor's relaxation. 

3.2.1. Abstraction of assignments. We focus on assignments (x%, . . . , xd) = T(xi,...,Xd) 
at control point i such that p o T is a quadratic function for every p £ P. Equation (3.1) 
translates in that case to (given that Vi-i defines the abstract value at control point i — 1, 
i.e. immediately before the assignment): 

{Ff(v))(p) = sup{poT(x) I q(x)-v i . 1 (q) < 0,Vg € P} . (3.7) 

We recognize the constrained optimization problem (3.2) and we use Lagrange duality as in 
the first step of Subsection 3.1 In our case, the Lagrange multipliers are some non- negative 
functions A from P to M.. We thus consider the function which we will call the relaxed 
function: 

(F?(v))(p) := inf suppoT(x) + T X(q) {v^q) - q(x)) . (3.8) 



AeJ-(P,M+) xG 



qeP 



To compute F^, we apply Shor's relaxation scheme of Subsection 3.1 and particularly the 
reformulation as the SDP problem (3.6), we get: 

(Fj l (v))(p)= Min n s.t. M(poT)+r ] N(-l) + y2X(q)(N(v^ 1 (q))-M(q)) < (3.9) 

\eJ r (P,R+) ^ p 



where M(poT), N(—l) and M(q) are the matrices defined in (3.5). 

We can treat the case where the map T has non-linear quadratic components. The 
condition p o T is quadratic for all p £ P implies, in this case, that the quadratic templates 
p should be linear forms which is equivalent to the fact that the matrices M(p) have the 
following form: 

i&f 
J& p 



(3.10) 
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and then the matrix ^^(—1) + X^eP A (#) (. v i-l{o)) — M(q)) in Equation (3.9) has no 



"quadratic counterpart" (of the form of Equation (3.10), we conclude that the SDP min 



imization problem is feasible in this case iff M(p o T) is negative semi-definite. Finally, 
the problem of Equation (3.7) is concave (since the forms q G P are linear) and conse- 



quently the Shor's relaxation scheme computes exactly the abstract semantic functional of 



Equation (3.7). 



x = [0,10]; 

y = i; [i] 

xn = — 3*x*x— y*y; 
yn = -y*y+x*x; 
x = xn ; 

y = yn [2]} 



Figure 6: A simple program written in polynomial arithmetic 



Example 3.3 (Non linear assignment). Let us take a simple program written in a polynomial 
arithmetic (assignments now involve polynomial expressions in the variables, rather than 
linear ones as in the previous example) . This program is described at Figure |6j We consider 
the set of templates P = {pi, P2}, with p\ : (x, y) 1— > x + y and P2 : (x, y) \-t x — y. We 
define the function T as follows: 

-3x 2 - y 2 
—y 2 + x 2 



(x,y) ^ T(x,y) 



which is the non- linear assignment of the program described at Figure pi By Equation (3.1 ), 

\2. 



the abstract semantics are, for w € J~(P, 

(FfH)(p) = (ll,9) 



(*i(«o)(p) 



sup (poT)(x,y) 
(x,y)e(i«i)* 



About the first component of the abstract semantics, the vector (11,9) means that 
(Ff (u>)) (p\) = 11 and (F^ (w)) (P2) = 9. This corresponds to the relations x + y < 11 and 
x — y < 9 at control point [1] . At control point [2] , the relation xn + yn = —2x 2 — 2y 2 holds 
between the variables at control point [1] and [2]. In this special case, the assignments are 
quadratic and the templates are linear, so the evaluation of the abstract semantic functional 
still reduces to a (non-convex) quadratic programming problem, which, by application of 
Shor relaxation, leads to the following relaxed functional (F^(w)^ (pi) at control point 2, 
for an element w G J-(P,M) 2 and for the template p\\ 

(^MKpi) = inf su p Yl x (q) w i(q) + ( x ,y) [ n 2 _°9 ) ( x ^y) T - Y] A (<?M^y) 

VgeP>0 i eF 



qeP 



By introducing the matrices 

(° 

Mfa) = 0.5 
\0.5 




M(p 2 
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M(pi o T) 




and using Equation (|3.9|), we get 



w))(pi)= Min 7] s.t. M(p 1 oT) + rjN(-l) + y^X(q)(N(wi(q))-M(q)) ^ . 
A(p)>o ^ 



Finally, using MatlatQ Yalmip [L04| and SeDuMi |Stu99j . we find (F^(w)) (jpi) = -3.2018e- 



09 ~ which means, at control point [2], that x < —y and the image of the second template 
by the relaxed function is (F^(w)) (p 2 ) = — 1.0398e — 09 ~ which means, at control point 
[2], x < y. The invariant found {(x, y) \ x < — y, x < y} is an unbounded set. We can refine 
this set from the invariant found by interval arithmetic to get a bounded set. 



3.2.2. Abstraction of simple tests. Now, we focus on a simple test and we write j the control 
point of the test. We assume here that a test is written as r(x%, . . . ,Xd) < where r is 
a quadratic function. We assume that the operation for the "then" branch has the form 
x = Tthen{ x ) and the operation for the "else" branch has the form x = T e i se (x) where Tthen, 
T e i se such that p o Tthen and p o T e i se are quadratic for all p G P. To enter into the "then" 
branch, the abstract values of control point j — 1 must satisfy the test condition of the "then" 
branch, so the abstraction of a test is Fj(X) = Tthen 

(Xj-i) n {x G R d | r(x) < 0} for the 
"then" branch. For the "else" branch, we have similarly Fj +2 (X) = T e / se (X,_i)n{x € R d \ 
r(x) > 0} which is equivalent to Fj +2 (X) = T eise (X,_i) n {x G R d \ -r[x) < 0}. However, 
we cannot use the Shor's relaxation scheme with strict inequalities and so we replace the set 
{x G ~R d | —r(x) < 0} by the set {x G M rf | — r(x) < 0} which is larger so we compute at least 
a safe over-approximation. When the function r is concave and the set {x G M rf | —r(x) < 0} 
is non-empty, the closure of the former set coincides with the latter set. For the "else" 
branch, the abstraction is, finally, Fj + 2{X) = T e / se (X/_i) n {x G M d | —r(x) < 0}. As we 
deal with arbitrary quadratic functions r, it is sufficient to show here how to deal with the 



equations at control point j and we simply write T instead of Tthen- By using Equation (3.1 ), 
we get, for v G F(P, M.) n , and p G P: 

(F)(v))(p) = (T n {x G R d I r(x) < 0})) (p) 

then, by a simple calculus: 

(F^(v)){p) = sup{poT(x) | q(x) < Vj-i(q) Vg G P, r{x) < 0}. (3.11) 
Using again the Lagrange duality, we get the relaxed problem: 

(*f(«))(p):= inf sup A( ? )D H (g) + poT(a) - rA(?)?(i) -^r(a). (3.12) 



4 Matlab is a registered trademark of the MathWorks,Inc. 
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Using the Shor's relaxation scheme described in Subsection 3.1 and in particular the SDP 



problem (3.6), we can rewrite (-F^(u)) (p) as the following SDP problem: 



Min n s.t. M(p oT) + r]N(—l) + X(q) (N(vi-i(q)) — M(q)) — uM(r) ^ (3.13) 
Ae^(P,M+) 



where M(poT), N(—l), M(q) and M(r) are the matrices defined in (3.5). 



x = [0 ,10] 
y = 1: 
u = 0; [1] 

if (y*y+x*x— 2<=0){ 
u=x ; 

x = 3-y*y; 
y = u-l;[2] 
} 

else { 

x=l-y*y ; 

y=3-y;[3] 

} 



Figure 7: A simple program written in polynomial arithmetic 



Example 3.4 (Quadratic assignment and quadratic tests). Let us take the simple program 
written in a polynomial arithmetic described at Figure [7j We introduce the quadratic 
function r : (x, y) h- > y 2 + x 2 — 2 which represents the test. We define the function T t h en as 
follows: 

(x, y) ^ T then (x, y) = f ^ _^ 
We also define the function T e / se by: 

l-y 2 

3-y 

These two functions represent the assignments of the branches then and else respectively. 



(x,y) H> T e i se (x,y) 



Similarly to Example 3.3, we consider the domain of intervals: the set of linear templates 



P = {x, —x, y, —y} where x:(x,y)t-tx ) y: (x, y) \-t y. The abstract semantics defined by 
Equation (j3TJ) are, for w G J"(P,IR) 3 : 

(Fl(w)){p) = {x(x,y) < 10, -x(x,y) < 0,y(x,y) < 1, -y(x,y) < -l} 1 " 

(F$(w)){p)= sup p{T then (x,y)) 

r(x,y)<0 

{ F h W ))(p)= SU P P( T else{x,y)) 
(x,y)£w* 
-r(x,y)<0 

At the first control point, for an element w € J~(P, M) 3 , the abstract semantics functional 
(i 7 ! (u>)) takes the value 10 for the template x, for the template —x, 1 for the template 
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y and -1 for the template —y. Again, the assignments are quadratic and the templates 
are linear, so the evaluation of the abstract semantic functional still reduces to a quadratic 
programming problem, which, by application of Shor relaxation, leads to the following 
relaxed functional (f£*(w))(x) at control point 2, for an element w G J-(P, M) 3 : 

{F?(w))(x)= inf sup^A^W^ + ^fU °)z-J2X{q)q(z) + 3-fir(z) 

By introducing the matrices 
M(x) 



r 


0.5 


0.5\ 




0.5 







, M(y) = j 


\0.5 











= -M(x) and M(-y) 


M(r) =10 1 I , M{xoT th , 



and using Equation (3.9), we get 





(F?(w))(x)= Min r, s.t. M(xoT fen )+7 ? AT(-l)-/iM(r)+y; A(g)(JV(«>i(g))-Af(g)) ^0 

A(p)>0 ^ 



We find using again Yalmip, SeDuMi and Matlab, (F^(w)^) (x) = 2. For the lower bound 
of the values of "3 — y 2 ", we find (-F^(w)) (— x) = —1. For the variable y at control point 
[2], we find {F^(w))(y) = -3.4698e - 09 ~ and (F^(w))(-y) = 1. We remark that 
(Ff-(w))(y) ~ which means that at control point [2], the values of "x — 1" are less or 
equal to 0. Indeed, Shor's relaxation scheme detects that the test is satisfied iff 1+x 2 — 2 < 
which is equivalent to x 2 < 1 and then the values of "x — 1" are bounded by 0. 



3.2.3. Properties of the relaxed semantics. We denote by A the set of coordinates of the 
abstract semantics functional which interpret an assignment. We denote by I the set of 
coordinates corresponding to tests (the symbol I stands for "intersection"). The set of 
coordinates U (for "union" ) corresponding to loops will be dealt with separately in Subsec- 
tion [3731 

Now, we are interested in the properties of the relaxed semantics which we introduced 



in the Subsection 3.2 First, we start by proving that the relaxed semantics is a safe over- 
approximation of the abstract semantics. It is deduced by the weak duality theorem: the 
relaxed semantics is a relaxation. 

Theorem 3.5. Let i be a coordinate in AUl. For all v G P{P, M) ra ; 

F*(v) < Ff-iv) 

The proof is given for the convenience of the reader. We only consider the case when 
i £ A (the proof can be easily adapted to the case of tests, i.e., 
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Proof. Let v G P(P, M) n and p £ P. If vf_ 1 is empty (this case includes Vi-\(p) = — oo for 
some p G P), (Ff(v))(p) = — oo for all p G P and the inequality holds. Now we suppose 
v*_i 7^ and we take x G Since A G J-(P,M+), E ge pA(g)(«i_i(g) - q(x)) > and 

then: 

V o T{x) < poT(x) + A(g)(uj_i(g) - g(z)) . 
<?eP 

We get: 

sup{po T(x) | < Vi-i(q), \/q G P} < sup po T(x) + X(q)(vi-\(q) - q(x)) . 

The left-hand side does not depend on A so we have [Ff(v))(p) < (F^(v)) (p). □ 

We observe that the calculus of the P-convex hull is a special case of assignment. Indeed, 
to compute a P-convex hull, is equivalent to compute the abstract effect of the assignment 
x = x. Consequently, we can apply the Shor's relaxation scheme to over-approximate the 
computation of the P-convex hull when the templates are quadratic. 

Corollary 3.6. Let w be in P(P, M) and p in P we have: 

(vex P (u/))(p) < inf 77 s.t. M(p) + r?iV(-l) + A(g) (N(w(q)) - M(q)) ^ 
\eJ r (P,R+) qep 



< w(p) . 



Proof. The first inequality is a special case of Theorem 3.5 Let us show the second inequal- 
ity. Let uu G P(P, M) and p G P. We know that the following equality holds: 

inf r] s.t. M(p) +rjN(-l) + ^ X(q)(N(w(q)) - M(q)) * 



= , jnf . SU P X] - V X(q)q(x) 

Let us choose X(q) = 1 if q = p and otherwise. Then, sup^g^d EoeP ^(l)' w ('l) + ^( 3; ) ~~ 
EoeP ^( ( l)l( x ) = 10 (p)- Since we have to take the infimum over A so the second inequality 
holds. □ 

The monotonicity of the relaxed semantics will be useful for the construction of Kleene 
iteration. Indeed, we will define in the next section a simple Kleene iteration so we show 
that the relaxed semantics constructed for assignments and tests define a monotone map. 
We adopt the convention that X(p)w(p) = if w(p) = 00 and A(p) = for some p G P. 

Proposition 3.7. For i G AU I, the map v \— > Fj^(v) is monotone on the set P(P, M) n . 

Proof. We only give a proof for the tests so let i be in I. Let v, w be in P(P, IR) n . We have 
for all p G P, the following equality: 

( F i l ( v )) (?) = in f yZ A (<?H'-i (<?) + sup p o T(x) - V X(q)q{x) - fir(x) . 

Since A G P(P, the map -U h-» EgeP A(?)«j-i(?) is monotone and since sup^gj^po 
T(x) — E«eP A ( ( ?)?( x ) ~~ f'f(x) is a constant (independent of i>) the function v 1— )• F^(v) is 
an infimum of monotone maps and thus a monotone map. □ 
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Now we focus on a property which will be useful to construct a policy iteration for 
quadratic zones. We are interested in the case in which the relaxed semantics (F^(y))(p) 
is equal to the supremum in Equation (3.8) for some A and equal to the supremum in 



Equation (3.12) for some pair (A,/i) in the case of tests. A simple condition provides the 



desired result: if, for all i G A U I, Slater constraint qualification holds i.e. there exists 
x G M d such that q(x) < Vi-i(q), Vg E P (and for this particular x, r(x) < holds for a 
test r), then there exists some A (and a couple (A,/x) is the case of tests) which achieves 
the minimum in (3.8) (the minimum in ( |3.12 ) in the case of tests). Moreover the over- 
approximation we make is not in general that big; in some cases, Inequality 3.5 is even an 
equality. 



Proposition 3.8 (Selection Property). 
ji G M+ such that: 



We assume that there exist A G J-(P,M + ) and 



sup p o T(x) + -%) { v i-i (<?) - Q( x )) ~ V r ( x ) 

is finite. If the set: 

{x G R d | q(x) - Vi-i{q) < 0, Vg G P] n {x G R d | r(x) < 0} 
is non-empty, then there exist A* G J-(P, M+) and \i* G M+ such that: 

(*f («))(p) = sup po T(x) + ^ A*(<7)(^_i(g) - g(x)) - /i*r(x) 

Furthermore, ifpoT is a concave quadratic form and if for allq G P such thatvi-x(q) < +oo, 
q is a convex quadratic form, then: 

(F] l (v))(p) = {Fl(v))(p) . 

Proof. The proof is given in the appendix. O 

In the cas e of assignment, we can reformulate a selection property as a particular case 
of Proposition 3.8 since we can replace fi by and remove the set {x G M. d | r(x) < 0}. We 



do not give the details for that. 



Remark 3.9. We can apply Proposition 3.8 in the case of intervals, zones and linear tem- 
plates: the functions p G P are all linear in these case s. When programs contain only linear 
expressions in assignments and tests, Proposition 3.8 implies that F^ = F^. Unfortunately, 
this seems to be the only simple case in which the relaxation is exact. In other words, the 
linear templates of Sankaranarayanan et al. are optimal in some sense. Indeed assume there 
is only one quadratic template p and consider the assignment y = T(x); T being affine. To 
evaluate the abstract semantic, we have to solve an optimisation problem of the form: 

Max p(y) 
s.t. y = T{x) 
p(x) < a 

Due to the constraint p{x) < a, the function p must be convex for the feasible set to be 
convex. However the same function also appears as the objective function and thus, p must 
be concave for the optimisation problem to be tractable by convex programming methods. 
This is possible only if p is affine. 
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The introduction of the duality provides a reformulation of Equation (3.8) and Equa- 



tion (3.12): the relaxed semantics can be rewritten as the infimum of affine functions of 
the function (when Vi-j(p) G K for all p £ P). Once again, this reformulation will be 
useful to construct a policy iteration since it allows to solve at each step a linear program. 
So, let us fix A G J r (P,M. + ) and observe that the sum Ylq^p ^(l) v i-i( < l) does not depend on 



the variable x in Equation (3.8). 



Let v be in F(P, R) n . We shall need the following notation. 
For i G A, we now define, for A G P(P,M.+), F*(v) by: 



{F l x (v))(p):=Tx( q )v^ 1 ( q ) + V^(p) . (3.14) 



qeP 



where V^{p) := sup po T(x) — X(q)q(x) (3.15) 

• For i G I, we define, for A G J"(P,M+) and y, G R+, ff'^t;) by: 

(i^(«))(p) := K<l)Vi-x{q) + V^(p) (3.16) 
qeP 

where V^'^ip) := sup p o T(x) — X(q)q(x) — yr(x) . (3-17) 
xm d q&P 

The relaxed functional can now be readily rewritten as follows. 
Lemma 3.10. For i G A and j G I: 

W(«))(p) = A6J igf R+) W(«))(p), (If («))(?») = Ae ; ( nf R+) ■ 



We remark that and ^ A,/i are the value of an unconstrained quadratic maximization 
problem. So, the functions V x and V^'^ 1 can be determined algebraically. Moreover, V^(p) 
and V X,fl (p) can take the value +oo if the matrices associated to the quadratic functions 
x i y p o T{x) — Ylq&p Mo)q( x ) an d x ^ po T(x) — J2 q eP ^{o)q( x ) ~~ fJ>r(x) are not negative 
semi-definite. Furthermore, the latter matrices depend on A G J~(P, M+) and on a couple 
(A,/i) G J"(P,K+) x K*. So, to ensure the finiteness of the value, it suffices to choose A (or 
a couple (A, y) in the case of tests) such that the corresponding matrix is negative semi- 
definite. We denote by A', the Moore-Penrose general inverse of A, which can be defined as 
lim^o A T (AA T + eld)^ 1 . The following proposition shows how to evaluate the functions 
V x and V i . We only consider the evaluation of V i since the evaluation of the former 
function can be viewed as a special case of the evaluation of the latter. We recall that all 
the assignments T are such that poT is a quadratic function for all p G P, and all functions 
r arising in tests are quadratic. Moreover, we recall that we write a quadratic function g as 
x i y x T A a x + b T n x + c . 



18 



A. ADJE, S. GAUBERT, AND E. GOUBAULT 



Proposition 3.11. Let i be inl and let an assignment T such that, for all q G P, qoT is 

a quadratic function. Let p be in P and let (A,/u) be a couple in J-(P, M+) x M + , we write: 



A p (X,fi) 


= Ap T - 


- 


— jiA r 






qeP 




£ p (a,m) 


= bpoT - 




\ib r 






qeP 




C p (\,n) 


= C po T - 




flC r 






q&P 





If(X,fi) G T(P,R+) x R + satisfies A(X, fi) ^ and B P (X, fi) G lm(A(X,n)) then: 

V^(P) = -\b p (X, y) T A p (X, fi)'Bp{\ M ) + C P (X, /i). 

Otherwise V i X,fl (p) = +oo. 

The proof is classical but it is provided for the convenience of the reader. 
Proof. Let (A, fj.) G P{P, M+) x M + and p G P. We can rewrite 

Vi'^ip) = SU P x T Ap(\, fi)x + B p (X,fj,) T x + C p (X,fi). 

We assume that A(X, y) is not a semi-definite matrix, there exists y G M rf and 7 G M, such 
that r y 2 y T A(X,/i)y > and jB p (X, y) T y > 0, now taking t > 0, we get t r y 2 y T A(X, \i)y > 
and t r jB p (X, y) T y > 0, this leads to T-^ > ty T A p (X,y)y + t / yB p (X, y) T y + C P (X, y) and 
we conclude that '^(p) = +00. 

We assume that A(X, y) < 0, the function x 1— > x T A p (X, fi)x + B P (X, y) T x + C P (A, //) is 
concave and differentiable then its maximum is achieved at every zero of its derivative. A 
maximizer x satisfies 2A P (X, y)x + B P (X, fi) = 0. If B P (X, y) £ lm(A(X,y)), the equation 
2A P (X, fi)x + B p (X, //) = does not have a solution. Taking a non-zero vector y in the kernel 
of A(X,fi) such that B p (X,fi) T y > 0, we get for all 7 > that V^(p) > jB p {X,fi) T y + 
Cp(A, fi). We conclude that V i X ' fl (p) = +00. Suppose that B P (X, y) G Im(A(X, fi)), then 

x = —-Ap(X, fi)'B p (X, y) + z where z belongs to the kernel of A P (X, y). Finally, we conclude 
that: 

V^( P ) = -\b p (X, y) T A p (X, ri'BpiX, n) + C P (X, n) 

since it suffices to compute x T A P {X, fj,)x + B p (X, y) T x + Cp(X, y) to find the value of V X,>1 (p). 

□ 



3.3. Abstraction of loops. The only point that we did not address yet is how to interpret 
the semantics equation at a control point i in which we collect the values of the variables 
before the entry in the body of the loop, at control point i — 1, with the values of the 
variables at the end of the body of the loop, at control point j: Fi(X) = Aj_i U Xj. By 
using Equation (3.1), for v G P(P,~S.) n , Ff(v) = (v*_i\Jvj)^ . As for zones, we notice that the 
union of two such P-convex functions Vi-\ and vj is directly given by taking their maximum 
on each element of the basis of quadratic functions P. Nevertheless, during the fixpoint 
iteration (as in SectionEl) the functions Vi-\ and Vj are not necessarily P-convex. Moreover, 
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if we take the abstract semantics FHv), we do not have an infimum of linear forms (or at 
least a maximum of linear forms) on the abstract values Vi-\ and Vj, a formulation that we 
need. Finally, we relaxed the abstract semantics fHv), using Remark 2.6 by the supremum 
itself and F^(v) = sup(uj_i, Vj). By this reduction, the map v i— > F^iv) is monotone on 
J~(P, M) n . Recall that U denotes the set of coordinates such that the concrete semantics is 
a meet operation. For i G U, the monotonicity of the map v i— > F^{v) follows trivially from 



the previous observations. Combining this with Proposition |3.7| above, we eventually get: 
Proposition 3.12. The map v i— > F^(v) is monotone on F(P,M.) n . 



To sum up, we conclude from Theorem 3.5 that we can compute over-approximations 



of (3.7) and (3.11) by solving a SDP problem. 

Theorem 3.13. In the case of quadratic templates, for a program with an affine arithmetics, 
the relaxed functional F^ can be evaluated using Shor's semi-definite relaxation and provides 
a sound over- approximation of the abstract functional FK 



4. Solving the semantic equation 

4.1. Fixpoint equations in quadratic zones. We recall that P is a finite set of quadratic 
templates. The map F is a monotone map which interprets a program with d variables and 
n labels in (V(M. d )) n . We recall that v* denotes the vector of sets ((ui)*, ,(v n )*) and 
F\v) = (F(u*))t i.e. Vi, f{(v) = (F i (w*))t and F n is the map, the components of which 
are the relaxed functions of F*. As usual in abstract interpretation, we are interested in 
solving the least fixpoint equation: 

M{v £ Vex P (P ^ W) n | < v} (4.1) 

Nevertheless, the function F^ is not easily computable (since the templates p are polyno- 
mials, the epigraph of F" can be checked to be a semi- algebraic set, but this of course does 
not lead to scalable algorithms). Hence, we solve instead the following fixpoint equation in 
F(P,R) n : 

mf{v e F(P,W) n | F n (v) < v} (4.2) 

and sometimes, we will stop our analysis at some vectors v such that F^(v) < v. 

We next describe and compare two ways of computing (or approximating) the smallest 



fixpoint of the semantic equation: Kleene iteration in Section |4.2[ and policy iteration in 
Section I 



4.2. Kleene iteration. We note by _L the smallest element of Vexp(P i— > M) n i.e. for 
all i = 1, • • • ,n and for all p E P, _Lj (p) = — oo. The Kleene iteration sequence in 
Vexp(P i— > M) n is thus as follows: 

(1) v° =!_ 

(2) for k > 0, v k+l = vexp oF n {v k ) 
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Since v i— > vexp oF^- is monotone over Vexp(P i— > M) n , this sequence is non-decreasing, and 
its limit is a candidate to be the smallest fixpoint of the functional vexp oF n . Unfortu- 
nately, it cannot be argued that this limit is always the smallest fixpoint without further 
assumptions. Indeed, F^, which is essentially defined as an infimum of affine functions, is 
automatically upper semi-continuous, whereas Scott continuity, i.e., lower semi-continuity 
in the present setting, would be required to show that the function F^ commutes with the 
supremum of increasing sequences. However, it can be checked that the map vexp oF K is 
concave (as the composition of a concave non-decreasing function, and of a concave func- 
tion), and it is known that a concave function with values in M is continuous on any open 
set on which it is finite |Roc96l Th. 10.1]. Hence, if the supremum of the sequence produced 
by the Kleene iteration belongs to such an open set, it is guaranteed to be the smallest fixed 
point. A more detailed theoretical analysis of the Kleene iteration, in the present setting 
of non-linear templates, appears to raise interesting technical convex analysis issues, which 
are beyond the scope of this paper. 

Kleene iteration has the inconvenience that the values v k which are obtained at a given 
iteration k (before convergence) do not provide a safe invariant. We shall see that policy 
iteration does not have this inconvenient: even if it is stopped at an intermediate step, it 
does provide a safe invariant. Moreover, the convergence of the Kleene iteration can be 
very slow, so it needs to be coupled with an acceleration technique which provides over- 
approximations. In our implementation, after a given number of iterations, and during a 
few iterations, we round bounds outwards with a decreasing precision (akin to the widening 
used in [GPBG08J). Note also that the P -convex hull cannot be computed exactly, so we 
over-approximate it using Shor relaxation. This yields an approximation of the sequence 
{v k )k>0i m which the approximated vectors v k do not belong necessarily to Vexp(P i— > R) n 
but only to F(P,W) n . 



4.3. Policy iteration algorithm. 

4.3.1. Selection property and policy iteration algorithm. A policy iteration algorithm can be 
used to solve a fixpoint equation for a monotone function written as an infimum of a family 
of simpler monotone functions, obtained by selecting policies, see |CGG + 05^ IGGTZ07] for 



more background. The idea is to solve a sequence of fixpoint problems involving simpler 
functions. 

In the present setting, we look for a representation of the relaxed function 

F n = inf F w (4.3) 
vren 

where the infimum is taken over a set n whose elements tt are called policies, and where 
each function F n is required to be monotone. The correctness of the algorithm relies on 
a selection property, meaning in the present setting that for each argument (i,v,p) of the 
function F 1 ^, there must exist a policy tt such that (Fj^(v)^j (p) = The idea of 

the algorithm is to start from a policy tt, compute the smallest fixpoint v of F w , evaluate 
F^ at point v, and, if v ^ F^{v), determine the new policy using the selection property at 
point v. 



Let us now identify the policies. Lemma 3.10 shows that for each template p, each 
coordinate Fj*~ corresponding to an assignment i 6 A can be written as the infimum of a 
family of affine functions v i— > F^(v), the infimum being taken over the set of Lagrange 
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multipliers A. The same lemma provides a representation of the same nature when the 
coordinate i £ I corresponds to a test, with now a couple of Lagrange multipliers (A,/u). 
Choosing a policy it consists in selecting, for each i G A (resp. j G I) and p G P, a Lagrange 
multiplier A (resp. a pair of Lagrange multipliers A,/u). We denote by 7Tj(p) (resp. 7Tj(p)) 
the value of A (resp. (A,/x)) chosen by the policy 7r. 

Then, the map -F 71 " in (4.3) is obtained by replacing F^ by the affine functions appearing 
in Lemma |3.10[ for i G A U I. For coordinates corresponding to loops, i.e., i G U, we take 
F? = (the choice of policy is tri vial) since the infimum operation does not appear in 
the expression of F n (see Subsection 3.3). 

Proposition |3.8| shows that the selection property is valid under a Slater constraint 
qualification condition. We thus introduce P5(P, M) n , the set of elements of P(P, M) which 
satisfy the Slater condition when the component Fi of F corresponds to an assignment or 
a test. More concretely: v G J-S(P,M) n , if, for all i G A the set: 

{x G M. d | q(x) < Ui_i(g), V<? G P} 

and, for i G I and a test r, the set: 

{x eR d \ q(x) < Vi 



V<? G P} n {x G M d I r(s) < 0} 



are non-empty. 



Algorithm 1 Policy Iteration in Quadratic Templates 



Choose 7T° G n, k = 0. 

Compute V 77 = {V 77 (q)} q ^p and defin e the associated f uncti on F n by choosing A and 
// according to policy 7r fc in Proposition 



3.11 



and Lemma 



Compute the smallest fixpoint v k in P(P,M)™ of F 1 " 
Compute w k = vexp(t> fc ). 

If w k G J-S(P, M) n continue otherwise return w k . 



3.10 



6 Evaluate F n (w k ), if F n (w k ) = w k 
F nk+1 (w k ). Increment k and go to 2. 



return w k otherwise take 7r fc+1 s.t. P^( 



Thi s lead s to Algorithm 111 For the third step of Algorithm [Tl since P is finite and using 

I — - — i 

F 71 " is monotone and affine P(P, M) n , we compute the smallest fixpoint of F 7 * 



3.10 



Lemma 

by solving the following linear program see [GGTZ07, Section 4]: 

n 

min E 1 («)) (9) — v k(<l)i V/c = 1, • • • , n, Vg G P (4.4) 

Remark 4.1. As in the case of the earlier policy iteration algorithms in static analy- 
sis lCGG+05l[GGTZ07j . an important issue is the choice of the initial policy, which may in- 
fluence the quality of the invariant which is eventually determined. In [ CGG + 05l IGGTZ07] . 
the initial policy was selected by assuming that the infimum is the expression of the 
functional is attained by terms corresponding to guard conditions, see specially § 4.2 
in [GGTZ07] . The same principle can be used here. Another method to choose an ini- 
tial policy is to run a few Kleene iterations, in combination with an acceleration technique. 
This leads to a postfixpoint v of P 7 ^, and we select as the initial policy any policy attaining 
the infimum when evaluating P^v) (i.e., choose for 7Ti(p) or TTj(p) any Lagrange multiplier 



A or pair of Lagrange multipliers A,/i attaining the infimum in Lemma 3.10) 
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Remark 4.2. To ensure the feasibility of the solution of (4.4) computed by the LP solver, 
we replace, when possible, the constraint set by F 71 ' (v) + e < v, where e is a small constant 
(typically of the order of several ulp(v), where ulp(v), which stands for "unit of least preci- 
sion" , is the minimum over the coordinates i of the differences between the nearest floating 
points around Vj). 

To obtain safe bounds even though we run our algorithm on machines which uses finite- 
precision arithmetic, we should use a guaranteed LP solver (e.g. LURUPA see |Kei05| ) to 
check that the solution obtained verifies F n (v) < v. 

In the fourth step of Algorithm[TJ the operation of closure is, in practice, the relaxation 
of the P-convex hull computed by a SDP solver (see Corollary 3.6). The same corollary 
shows that the SDP relaxation of the P-convex hull of some w G J-(P,M), is still smaller 
than w and this result ensures a gain of precision. 

We can only prove that policy iteration on quadratic templates converges (maybe in 
infinite time) towards a postfixpoint of our abstract functional and that under some tech- 
nical conditions, it converges towards a fixpoint. One interest in policy iteration for static 
analysis is that we can always terminate the iteration after a finite time, and end up with 
a postfixpoint. 

Theorem 4.3. The following assertions hold: 

(1) F n (v l ) £ v l F n (v l ) < v l ; 

(2) The sequence v l computed by Algorithm^ is strictly decreasing; 

(3) The limit v°° of the sequence v l is a postfixpoint: F^(v°°) < v°° . 

Proof. (1). Let I G N. We assume that I > and F n {v l ) / v l , there exists ir l such that, 
F %l (v l ) = v l and since F n = inf F w , we get F n (v l ) < F n \v l ) = v l and from F n (v l ) / v\ 
we conclude that F^(v l ) < v l . 

(2) . We prove the second point by induction on I G N. We suppose that I = if 
F^(v°) ^ v° and that v° G FS(P, M) n otherwise the algorithm stops. There exists it 1 
such that F n (v°) = F w (v°) < v° by the point 1. Moreover, v 1 is the smallest element of 
{v G F(P, M) n | F n (v) < v} thus v 1 < v° and since F n (v°) < v° we conclude that v l < v°. 
The same argument holds if v l < v 1 ^ 1 and v l G FS(P, M) n . 

(3) . Finally, we deduce from v°° < v that F^^v 00 ) < F^{v ) < v\ Taking the infimum 
over /, we get F n (v°°) < v°°. □ 

Remark 4.4. It is desirable to choose (if possible) the initial policy ir° so that the set: 

{v£ F(P,R) n \ (Ff'(v))(q)<v i (q), Vi = 1, • • • , n, Vg G P} 
is non-empty. Indeed, the non-emptyness of this set e nsure s that the coordinates of the first 



vector v° are not equal to +oo and then, by Theorem 4.3, all the terms w k of the sequence 
generated by the policy iteration have coordinates which are not equal to +oo. Then, 
the policy iteration algorithm at any step k and at any breakpoint i returns non-trivial 
invariants of the form q{x) < a (with a := wf(q) finite). 

Remark 4.5. The policy iteration algorithm developed in |CGG + 05[ IGGTZ07] can be re- 
covered as a special case of Algorithm [TJ when applied to a domain of linear templates or to 
the domains of zones or intervals, for a program containing only linear expressions in assign- 
ments and tests. Indeed, the main addition in the present algorithm is the presence of the 



relaxation, and the latter turns out to be exact in these special cases, see Proposition 3.8 
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X = [0 ,1]; 

v = [0,1]; [1] 

h = 0.01; 

while (true) { [2] 
u = v; 

v = v*(l — h)— h*x; 
x = x+h*u; [3] } 



Ff(w){p) = {x(x,v) < 1, v(x,v) < 1, L(x,v) < 7} 
F k w )(p) = (sup{wJi(p),^(p)}) t 
F^(w)(p) = sup p(T(x,v)) 

Figure 8: Implementation of the harmonic oscillator and its semantics in J-(P, M) 3 

4.4. Max strategy iteration. Gawlitza and Seidl [GSlOj developed an alternative iter- 
ation to compute the least fixpoint of the relaxed semantics for the quadratic templates. 
The relaxed semantics which they use are constructed from the dual program of Shor's 



relaxation SDP problem (3.6). Our relaxed semantics coincide with their relaxed semantics 
when technical conditions (which are often satisfied) hold. They obtain a map whose coor- 
dinates are the maximum of a finite number of concave functions. Their approach consists 
in solving the least fixpoint equation from below and they initialize their iteration by the 
function which is identically equal to — oo. At each step of their algorithm, they select a 
function which achieves the maximum. Since they compute the least fixpoint from below, 
the least fixpoint is returned. This implies that their iteration must be run until a fixpoint 
is reached whereas our approach allows to stop the iteration at each step of the algorithm 



to provide a valid invariant. A survey GSA + j recapitulates the two approaches 



4.5. A detailed calculation on the running example. Now we give details on the 
harmonic oscillator of Example [TJ The program of this example which is given at Figure [8] 
implements an Euler explicit scheme with a small step h = 0.01, that is, which simulates 
the linear system (x, v) T = T(x, v) T with 

1 h 
-h 1-h; 

We want to use the information of a Lyapunov function L of the linear system T to compute 
bounds on the values taken by the variables x and v of the simulation: the function (x, v) :\— )■ 
(x, v)L(x, v) T furnishes a Lyapunov function with 

'2 r 



T 



L 



1 3 



We also use the quadratic functions (x, v ) h-> x 2 and (x, v) t— > v 2 which corresponds to 
interval constraints. We introduce the set of templates P = {x, v, L} and below the program 
it is described the semantic equations for all the three control points. Now we are going to 
focus on the third coordinate of (F^(v)^(p). Let us consider, for example, p = x, we get: 
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inf sup X(q)w2(q) + (ay 



1 - A(x) h/2 
h/2 h 2 - X(v] 



X(L)L)(x,vf (4.5) 



By introducing the following symmetric matrices, we can rewrite (4.5) as Equation (4.6): 

'o o o\ /o o o N 

I and M(x o T) 



M(x) 



10 
,0 0, 



M(v) 



.0 1. 




h/2 h 1 



(F^(w))(x)= Min rj s.t. M(xoT)+rjN(-l)+ ^ X(q) (N(w 2 (q)) - M(q)) < 

^^R^ q=x,v,L 

(4-6) 

To initialize Algorithm 1 1 we choose a policy 7r°. For the third coordinate of F^, we have 

to choose a policy 713 such that V 3 3 (p) is finite for every p = x,v,L. We can start, for 
example, by: 

vr°(x) = (0,0,1), 4(v) = (0,0,1), tt§(L) = (0,0,1) . 
This consists, for p = x, in taking A(x) = X(v) = and A(L) = 1 in (4.5). By Proposition 
l3~TTl we find: 

, . / -1 /i/2-T 

■ r -' ] {h;2 1 /. 2 -3 l, ' r - r) 



V^ 3 (x) = sup (av 



3 (v) = sup ( ^ )eR2 (x,i;) L^^vj , 



/, 2 



/i(l - h) - 1 
(1 -/i) 2 -3 



(x, v) 



T 







Vf 3 (£) = sup {XtV)m 2(x, v)(T T LT - L)(x,v) T = 

The solution of the maximization problems are zero since all the three matrices are negative 
definite (i.e. a matrix B is negative definite iff x*^4x < for all 1 / 0). The third matrix 
T T LT — L is negative definite since L satisfy the Lyapunov condition for the discrete linear 
system (x,v) = T(x,v). To compute the least fixpoint of F w , we solve the following linear 
program (see (4.4)): 

3 

h{L)<p 3 (x), i3 2 (L)<f3 3 (v), fS 2 (L)</S 3 (L) 
fefe)</3 2 fe), Pz{v)<l3 2 {v), Pz{L)<f3 2 (L) 

l<fa(x), l</3 2 fe), 7</3 2 (i) 

l<Afe), l<Pi(v), 7</3i(L) 

Using solver Linprog, we find: 

u?(x) = 1.0000 u Q 2 {x) = 7.0000 u§(s) = 7.0000 
u{(v) = 1.0000 u§(u) = 7.0000 uf(v) = 7.0000 
u?(L) = 7.0000 u§(L) = 7.0000 u§(X) = 7.0000 
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{x 2 < 4.2000, v 2 < 2.8000, 2a; 2 + 3v 2 + 2xv < 7} 

{x 2 < 4.1723, v 2 < 2.7724, 2a; 2 + 3v 2 + 2xv < 7} 

{a; 2 < 3.9386, v 2 < 2.5940, 2x 2 + 3v 2 + 2xv < 7} 




{a: 2 < 3.6079, v 2 < 2.3876, 2a; 2 + 3v 2 + 2xv < 7} 

{a- 2 < 3.5051, v 2 < 2.3353, 2a- 2 + 3v 2 + 2xv < 7} 

{x 2 < 3.5000, v 2 < 2.3333, 2x 2 + 3v 2 + 2xv < 7} 

Figure 9: Successive templates along policy iteration, at control point 2, for the harmonic 
oscillator. 

The approximation of the closure of u°, vexp(ii°) is given by a Matlab implementation, 
using Yalmip and SeDuMi returns the vector w\: 

wl(x) = 1.0000 w° 2 {x) = 4.2000 w%(x) = 4.2000 
wl(v) = 1.0000 w{(v) = 2.8000 wf(v) = 2.8000 
wl(L) = 7.0000 w{(L) = 7.0000 wf(L) = 7.0000 

Using again Yalmip with the solver SeDuMi, the vector w is not a fixpoint of F^, so we get 
the new following policy: 

ttI(x) = (0,0,0.596), irl(v) = (0,0,0.3961), t^(L) = (0,0,0.9946) . 

Finally, after 5 iterations we find that the invariant of the loop i.e. at control point 2 is 
the set: 

{x 2 < 3.5000, v 2 < 2.3333, 2a; 2 + 3v 2 + 2xv < 7} . 
We draw w 2 at each iteration of Algorithm [T] in Figure [9j 

This method is to be compared with the classical Kleene iteration with widening. On 
this example, we find without widening x 2 < 3.5000, v 2 < 2.3333 and 2x 2 + 3v 2 + 2xv < 7 



in 1188 iterations whereas with the acceleration technique described Subsection 4.2 we find 
x 2 < 6.0000, v 2 < 4.0000 and 2x 2 + 3t; 2 + 2xv < 10 in 15 iterations. 

5. Benchmarks 

We implemented an analyzer for the quadratic template domain we presented, written in 
Matlab version 7.8(R2009a). This analyzer takes a text file in argument. This text file 
corresponds to the abstract equation v = F$(v) where F$ is defined by Equation (2.4). The 
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Figure 10: Benchmarks 



quadratic template can be loaded from a dat file by the analyzer. The afhne maps are 
treated in the same manner. 

In this analyzer, we can choose to use the Kleene iteration method or policy iteration. 
For the Kleene iteration method, the user gives as an argument a maximal number of 
iteration and if the acceleration method has to be applied. The acceleration method start 
from the iteration n + 1 (n denotes the number of lines of the code) and ends eleven iteration 
after. For the policy iteration method, the user gives the dat file defining the initial policy 
or chooses to make Kleene iterations before determining the initial policy. 

For the policy iteration, the user gives also as argument a maximal number of iteration 
and the policy iteration runs until a fixpoint is reached or Slater constraint qualification is 
no longer satisfied or if the maximal number of iteration is achieved and so a postfixpoint 
can be returned by the policy iteration algorithm. Similarly, the Kleene iteration with 
acceleration provides a postfixpoint after acceleration and widening to top, if the iteration 
does not converge after a given number of iterations. The analyzer writes, in a text file, 
information about time, quality of the invariants found and number of iterations. 

For the benchmarks, we used a single core of a laptop PC Intel(R) Duo CPU P8600 at 
2.4 Ghz with a main memory of 4Gb. We indicate in Table [TO] , the name of the program 
analyzed, the method used (policy iteration or Kleene iteration) for solving the fixpoint 
equation, the cardinality of the basis of quadratic templates used, the number of lines 
of C code the program has, the number of variables it manipulates, the number of loops. 
Then we indicate the number of iterations made, whether it reaches a fixpoint or (strictly) a 
postfixpoint. Finally, the last column concerns the time in seconds to compute the invariant. 
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x = [0 ,1]; 

Y = [0,1]; [1] 
while [2] (true) { 

x = (3/4)*x-(l/8)*y; 

y = x; [3] 

} 



Figure 11: The program Filter 

The file RotationlO is the problem of Example |3.2| in dimension 10. By the fixpoint 
computation, we prove automatically that the unit sphere in dimension 10 is invariant by 
rotation. Both Kleene iteration and policy iteration find the unit sphere as invariant. 

The program Filter is an implementation of recursive linear filter of second order. The 
program is described Figure [TTj By policy iteration, we find the following set as loop 
invariant (at control point 2): 

{-0.5 < x < 1, -0.5 < v < 1, 3x 2 + v 2 < 4} . 

Although, the Kleene iteration with acceleration finds the following set: 

{-1.8257 < x < 1.8257, -3.1623 <v< 3.1623, 3x 2 + v 2 < 10} . 

The program Oscillator is the problem [TJ The invariant depicted Figure [TJ in Section [TJ is 
found by policy iteration whereas Kleene iteration after applying acceleration techniques 
from the iteration 4 to iteration 15 finds the less precise invariant {x 2 < 6.0000, v 2 < 
4, 2x 2 + 3v 2 + 2xv < 10}, in more time. 

In order to illustrate the scalability of the method, we considered a higher dimensional 
analogue of the example of Figure [TJ modelling N coupled harmonic oscillators, Xi + Xi + 
Xi + e^i<j<Ar x j = 0, for 1 < i < N. Hence, ./V variables Xi, Vi, and Wi now appear in the 
discretised scheme instead of x,v,w. We computed automatically a Lyapunov function in 
Matlab. We used as templates the one arising from the latter Lyapunov function, together 
with YliLi x i an d v h We took e = 0.5. We made tests successively for N = 2, 5, 10, 20 
(Oscillatorc2, Oscillatorc5, OscillatorclO and Oscillatorc20). As in the program Oscillator, 
we were interested in the loop invariant. For example, for N = 20, we found with policy 
iteration algorithm the following set: 

20 20 

{(x,v) E M 20 x M 20 | J^xf < 115.7169, ^ v 2 < 266.3048, (x, v)L(x, v) T < 350.0690} 

i=l i=l 

where L denotes positive semi-definite matrix associated to the Lyapunov function found 
automatically by Matlab. 

The Symplectic example implements a discretisation of x + cx + x = with c = by 
a symplectic method, considering specially the case in which c = (there is no damping). 
Then, the dynamical system has imaginary eigenvalues (its orbits are circles), reflecting 
the fact that energy is constant. However, the Euler scheme, which does not preserve 
this conservation law, diverges, so we use a symplectic discretisation scheme (preserving 
the symplectic form, see [HLW03J). This is an interesting, highly degenerate, numerical 
example from the point of view of static analysis, because there is no "stability margin", 
hence, methods not exploiting the Lyapunov function are likely to produce trivial invariants 
when c = 0. As in Oscillator, we start from a position x S [0, 1] and a speed v G [0, 1]. 
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tau = 0.1; 




x = [0 ,1]; 




v = [0,1]; [1] 




while [2] (true) 


{ 


x = (l-(tau/2) 


)*x+(tau -((tau ~3) /4))*v; 


v = — tau*x+(l- 


(tau/2))*v; [3] 


} 



Figure 12: An implementation of the symplectic method 



The discretisation of x + x = with the symplectic method and a step r = 0.1 gives us the 

3 

matrix T such that T\ y \ = 1 — |, = r — T2 j i = — r and T2 5 2 = 1 — |. We use the 
Lyapunov function L such that L(x,v) = (x,v)Q(x,v) with 

'1 



Q = 



o i-£ 



The symplectic method ensures that L(T(x,v)) = L(x,v). Our method takes advantage of 
this conservation law, since L is embedded as a template. The policy iteration returns the 
following as invariant set at the control point 2: 

{-1.41333 < x < 1.41333, -1.4151 < v < 1.4151, x 2 + 0.9975t; 2 < 1.9975} . 

Whereas the Kleene iteration returns: 

{-3.16624 < x < 3.16624, -3.16628 < v < 3.16628, x 2 + 0.9975w 2 < 10} . 

which is less precise. In particular, the Kleene algorithm misses the invariance of the 
Lyapunov function. 

SymplecticSeu is a symplectic method with a threshold on n = x. We iterate the 
Symplectic method while v > |, which gives the following code: 



x = [0,1]; 

v = [0,1]; 

tau = 0.1 [1] 

while [2] ((v> = l/2)) { 

x = (1 — tau /2)*x+(tau — (tau " 3) /4) * v ; 

v = — tau*x+(l — tau/2)*v ; [3] 

}; 



The policy iteration returns the following set which the invariant found at control point 2: 

{0 < x < 1.3654, < v < 1, x 2 + 0.9975w 2 < 1.9975} . 

and the policy iteration returns the following set at control point 3: 

{0.0499 < x < 1.3222, 0.5 < v < 0.9950, x 2 + 0.9975v 2 < 1.9975} . 

Although, the Kleene iteration with acceleration provides the following set which the in- 
variant found at control point 2: 

{0 < x < 3.1623, < v < 3.1662, x 2 + 0.9975w 2 < 10} . 

and the same set at the control point 3. 
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Remark 5.1. In the present benchmarks, the execution time of a single policy iteration step 
and of a single Kleene iteration step are comparable (the acceleration provided by policy 
iteration comes from the smaller number of iterations). Indeed, at each Kleene or Policy 
iteration, the bottleneck appears to be the evaluation of the relaxed functional F^, which 
requires to solve a family a Shor SDP relaxations. Policy iteration requires in addition to 
solve a family of linear programs (to compute the smallest fixpoint of the current policy). 
This is technically easier than solving the SDPs, and generally faster. However, each SDPs 
is typically local (involving a small number of variables), whereas the linear programs, which 
couple all the breakpoints and all the templates functions, may be of large size. Hence, one 
may construct (very large) instances in which solving the linear programs would become 
the bottleneck. 

6. Conclusion and future work 

We have presented in this paper a generalization of the linear templates of Sankaranarayanan 
et al. [SSM05j SCSM06J allowing one to deal with non- linear templates. We showed that 
in the case of quadratic templates, we can efficiently abstract the semantic functionals 
using Shor's relaxation, and compute the resulting fixpoint using policy iteration. Future 
work include the use of tighter relaxations for quadratic problems. In particular, sum of 
squares (SOS) relaxations (see for instance |Las07| and |Par03| ) would allow us to obtain 
more accurate safe over-approximation of the abstract semantic functional for arbitrary 
polynomial templates and for a general program arithmetics. Kleene iteration could be 
easily implemented in that way. However, the issue of coupling SOS relaxations with policy 
iteration appears to be more difficult. Note also that unlike Shor relaxation, SOS relaxations 
are subject to a "curse of dimensionality" . These issues will be examined further elsewhere. 
Another problem is to extend the minimality result of [AGG08 which is currently only 
available for the interval domain, to our template domain. We intend to study more in- 
depth the complexity issues raised by our general policy iteration algorithm. Finally, we 
note that a more detailed account of the present work has appeared in the Phd thesis of 
the first author |Adjll| . 
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Appendix 



Let / be a function from R d to MU {+00}. We say that the function / is proper iff there 
exists x £R d such that f(x) G R. 

We recall that / is convex iff the set epi(/) := {(x,a) £l d xl f(x) < a} is a convex 
set (i.e. for all t G [0, 1], for all x,y G epi(/), tx + (1 — t)y G epi(/)). Since the intersection 
of a family of convex sets is a convex set then the pointwise supremum of a family of convex 
functions is also a convex function. 

The function / is said to be lower semi-continuous iff the sets {x G M. d \ f(x) < a} are 
(topologically) closed for all a£l. Since the intersection of a family of closed sets is closed 
set then the pointwise supremum of a family of lower semi-continuous functions is also a 
lower semi-continuous function. 

The function / is level-bounded iff {x G M. d \ f(x) < a} are bounded set for all a£l. 

If lim f{x) = +00 then / is level-bounded. 

||x||-H-do 

Proposition |3.8| is deduced from a well-known result of convex optimization: 

Proposition 6.1. Let f : M. d \— > MU{+oo} be a convex, proper, lower semi- continuous and 
level bounded function. Then inf f(x) is finite and there exists x such that: 

f(x) = inf f(x) . 



A proof of Proposition 6.1 can be found in [AT03, Proposition 3.1.3]. 

Hr{x) if A G.F(P,] 



Proof of Proposition 3.8. We begin the proof by writing: 
sup poT(x) +^2\(q)(vi-i(q) - q(x)) 



qeP 



+00 



f ),AtGM+ 
otherwise 



To show Proposition 3.8 it suffices to show that the hypothesis of Proposition 6.1 



are 



verified: the function g is convex, lower semi-continuous, proper and level-bounded. 

The set P is finite so P(P,M.) is the finite dimensional vector space where \P\ is 
the cardinality of P. First, the function g is convex and lower-continuous as the pointwise 
supremum of convex continuous functions. 

The function g is also proper since g(\,(i) > —00, for all A G J-(P,M.+ ), for all [i G M+ 
and there exist A G T{P y M+) and fj, G M+ such that g(X, /j.) < +00. 

The function g is level-bounded since A G J-(P,M.+), [x G M+ and there exists some 
x G M. d such that, for all q G P, Vi-i(q) — q{x) > and r(x) < 0, we conclude that 
g(X, /i)>po T(x) + mm qeP (vi-i(q) - q(x)) J2 qe p •%) ~ r(x)(i then g(X, n) tends to +00 
as maxdlAj^, |/x|) tends to 00. 

The second part of Proposition |3.8| follows from the strong duality theorem for convex 
optimization problems, see e.g. [AT03, Proposition 5.3.2]. □ 



This work is licensed under the Creative Commons Attribution-NoDerivs License. To view 
a copy of this license, visit http://creativecommons.0rg/iicenses/by-nd/2.o/ or send a 
letter to Creative Commons, 171 Second St, Suite 300, San Francisco, CA 94105, USA, or 
Eisenacher Strasse 2, 10777 Berlin, Germany 



